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We study strongly correlated Hubbard systems extended to symmetric A^-component fermions. 
We focus on the intermediate-temperature regime between magnetic superexchange and interac- 
tion energy, which is relevant to current ultracold fermionic atom experiments. The A/'-component 
fermions are represented by slave particles, and, by using a diagrammatic technique based on the 
atomic limit, spectral functions are analytically obtained as a function of temperature, filling factor 
and the component number N . We also apply this analytical technique to the calculation of lattice 
modulation experiments. We compute the production rate of double occupancy induced by modu- 
lation of an optical lattice potential. Furthermore, we extend the analysis to take into account the 
trapping potential by use of the local density approximation. We find an excellent agreement with 
recent experiments on ""^^"^Yb atoms. 
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I. INTRODUCTION 

The realization of ultracold atoms in an optical lat- 
tice opens up the possibility to study in a controlled 
way strongly correlated quantum systems [l|, |2|. Such 
strongly correlated atoms are well described by the Hub- 
bard model. This model plays a central role for the study 
of the Mott insulator (MI) transition [s], high- Tc super- 
conductivity [4], and quantum magnetism. In addition, 
the high controllability of model parameters such as the 
interaction by a Feshbach resonance technique and ki- 
netic energy by changing the lattice depth allows us to 
capture such Hubbard physics in a broad range of pa- 
rameter regimes. 

The recent achievement of Fermi degeneracy of ultra- 
cold alkaline-earth- metal (-like) atoms such as ^^Ca, ^^Sr, 
and ^'^^Yb potentially provides a new class of strongly 
correlated matter. The structure of their nuclear spin de- 
grees of freedom allows the realization of high symmetry 
groups for the internal degree of freedom ("spin"). For 
instance, ^''^Yb atoms behave as a spin-5/2 fermion 
[7I] and ^^Sr as a spin-9/2 fermion [sHllj. In particular, 
provided all 5- wave scattering lengths are independent of 
the atomic spin states, the atom cloud as a many-body 
system obeys high symmetries. Thus the confinement 
of alkaline-earth- metal (-like) atoms in an optical lattice 
provides opportunities for the study of the SU(A^) sym- 
metric Hubbard model for spin-(A^ — l)/2 atoms. The 
experimental realization of such SU(A^) Hubbard mod- 
els has stronglystimulated the corresponding theoretical 

studies BliEIiMil. 

In condensed matter physics, the SU(A^) symmetric 
systems have been introduced as a purely theoretical ex- 
tension of the strongly correlated electron systems with 
SU(2) spin rotational symmetry, e.g., for quantum mag- 
netism liO-Iil or for the Hubbard model [25| in the 
context of high- Tc superconductivity [26-29|. As a the- 
oretical tool to solve such problems, the slave-particle 
technique, originally developed for the single impurity 



Anderson model |3Ql-l33ll, has been used and applied to 
the Hubbard model [34-38]. More recently, the slave- 
boson approach has also been used with success in the 
cold atom context (SOMiH. 

In this paper, we generalize the slave-boson calcula- 
tion scheme introduced in Ref. [41] to the symmetric 
A/'-component fermionic atom systems including SU(A/') 
symmetry. This technique has proven an effective way 
to compute the dynamics of strongly interacting systems 
at a filling of one or less than one particle per site in the 
spin-incoherent temperature regime for which the tem- 
perature is lower than the interaction energy, but larger 
than the magnetic superexchange one. This regime is di- 
rectly relevant to the current experiments on fermionic 
atoms. A diagrammatic approach based on the noncross- 
ing approximation (NCA) with the spin-incoherent as- 
sumption is used to estimate self-energies and compute 
the spectral functions as functions of temperature, chem- 
ical potential, and component number A^. This allows us 
in particular to compare the physics for different A's. 

We also use these techniques to compute the effect 
of lattice modulation spectroscopy which has been re- 
cently implemented in experiments. The lattice modu- 
lation technique has been originally applied to bosonic 
atom systems, in which the absorbed energy is measured 
as a function of the modulation frequency Ac- 
cording to the linear response formalism, the energy ab- 
sorption rate in such a weak perturbation regime gives 
access to the kinetic-energy correlation functions (isl - 
For fermionic atoms, an accurate measurement of 
the absorption energy is difficult, and a variant of the 
probe measuring the production rate of so-called dou- 
blons, which is the number of doubly occupied sites, in- 
duced by the lattice modulation has been proposed. [45| 
The doublon production rate (DPR) has been shown to 
be identical to the energy absorption rate both numer- 
ically and analytically [45-47]. This doublon measure- 
ment technique has been successfully implemented in a 
fermionic atom experiment [48]. This allowed more re- 
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cent experiments to successfully reach the linear response 
regime for which the DPR spectrum scales quadratically 
with the modulation amplitude. [49| A direct comparison 
with equilibrium theory is possible and has been success- 
fully done [39, 41]. So far the lattice modulation experi- 
ment has been done with ^^K [49] behaving as a spin- 1/2 
fermion and ^^^Yb [50] behaving as a spin-5/2 fermion. 

The paper is organized as follows. In Sec. |TT1 we intro- 
duce the Hubbard model for the A^-component fermions 
which includes SU(A') symmetry. The introduced Hamil- 
tonian is rewritten in a slave-particle representation. In 
Sec. lIIH we discuss the single particle properties based on 
the Hubbard Hamiltonian using this slave-particle repre- 
sentation. Then, using a diagrammatic approach starting 
from the atomic limit, self-consistent equations of dou- 
blon and holon self-energies are constructed and solved 
analytically. In Sec. lIVi using the spectral functions given 
in Sec. ITTTl we investigate the DPR spectrum produced 
by an amplitude modulation of the optical lattice po- 
tential. Then the analytic form of the DPR spectrum 
is given. In addition, we extend the calculation to the 
trapped case by using the local density approximation 
(LDA), and we also compare our results to the recent 
experiment with ^^^Yb atoms. Finally, our results are 
summarized in Sec. |Vl Some technical details on the for- 
mulation of the DPR spectrum are briefly reviewed in 
Appendix El 



II. MODEL 

A. SU(A)-symmetric Hubbard model 

In lattice systems with multicomponent particles, mul- 
tiparticle occupation states in addition to double occu- 
pancy can be defined in general. However, when the in- 
teraction between different components is strong, such 
multiple-occupation is at a higher energy state than dou- 
ble occupancy. Because we are interested in physics of 
doublon excitations at a filling of one or less than one 
particle per site, such higher occupation states are way 
above the main energy scale of interest. Thus, as an effec- 
tive model Hamiltonian, we can extend the Hubbard type 
two-body interaction to the A^-component case. Then in- 
teractions between different components, determined by 
the 5- wave scattering lengths, generally take different val- 
ues depending on the components: The interaction pa- 
rameter is written as Ua,a' where a and a' are the indices 
characterizing the internal state of the fermions. We con- 
sider the s pec ial case of a unique interaction parameter, 
Ua,a' = U [HI, [52); namely the coupling does not depend 
on the components. [53] Then, the interaction term has 
the same symmetry as the kinetic term, and the system 
symmetry turns out to be enlarged to SU(A/') symmetry. 

We consider the generalized A^-component fermionic 



Hubbard model, H = Hk + i^at with 



N 



^=1 (iJ) 



N 



(1) 



a=l j 



where Cj^cr is the annihilation operator of a fermion with 
the internal component a at a site j, and nj^a- is the num- 
ber operator. The parameters J and [/, respectively, de- 
note the nearest-neighbor hopping energy and the on-site 
interaction between components a and (T^ Throughout 
this paper, we consider only the repulsive case U > 0. 

In the considered regime of chemical potentials, 
particle-hole symmetry always disappears except for the 
N = 2 case at /i = [7/2. This is because for A" > 2 there 
are A" (A" — l)/2 doublon states while the hole state is 
unique. 



B. slave particle representation 

The A^-component fermions have a larger Hilbert space 
than that of the two-component case: While an empty 
site (holon) is unique, there exist multiple-occupation 
states (three- and four-fold occupation, and so on) in 
addition to A" (A" — l)/2-species doubly occupied states 
(doublons) and A^-species single occupied spin states 
(spinous). The multiple-occupation states are energet- 
ically out of shell, since those states cost an energy 
higher than 2[/, which is over the energetic cutoff in 
our model Hamiltonian. Thus in our case we describe 
the single-site state by a holon A'-species spinon 
|cr) (cr = 1,2, • • • , A'), and A'(A' - l)/2-species doublons 
\da,a') {o' ^ cr' and a^a^ = 1, • • • , A^). In terms of dou- 
blon states, the antisymmetrization condition \dcr,a') = 
— \da',a) is imposed. We hereafter suppose that the spin 
indices in \da,a') are ordered such that a > c^^ The 
single-site original fermionic operators, Ca- and c^, are 
given by 

cr-l N 

c. = \h){a\+Y,W'){d.,.'\- E (2) 

cr-1 N 

cl = \a){h\+J2K,,,){a'\- E K..')(<t'|. (3) 

a' = l a'=a-\-l 

Then the representation no longer gives back the anti- 
commutation relation {ccr,cj.} = 6a,a'^ but it should be 
approximately correct as long as U is much larger than 
the particle hopping J and temperature, and the filling 
is less than unity, which means that all the excitations 
leave the system in the proper subpart of the Hilbert 
space. We introduce the creation and annihilation oper- 
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ators for a holon, spinons, and doublons as 



method. The Hamiltonian ([T]) is represented as 



\a) {h\ = hlK 
\h){cj\=h^h,, 



and the new vacuum state |0) is defined as |0) = 
h\{))=d,,> |0) = 0. 

It is easy to extend the above single-site argument 
to multi-site problems: All operators become site de- 
pendent. In order to recover the anticommutation re- 
lations between the original fermions at different sites, 
{ci,cr, ^\ct'} — ^ ^ j)^ we assume the following commu- 
tation and anticommutation relations: {hi^h^j} = Sij^ 

{di^aa'.d'^j^^r}'} = ^i,i^a,r7^a' ,77' , and ^],^/] = ^i,j^a,a'' 

Furthermore, by imposing the following constraint we 
project onto the physical subspace the Hilbert space en- 
larged by introducing the holon, spinons, and doublons: 



AT 



N 



.cr=l 



H.c 



E 



TrC^3Cr2jt J. 



Cl>C2>cr3 



H.C. 



N 



3 



cr>cr' 



(8) 



(9) 



where the local potentials for the slave particles have 
been introduced as 



b 
3 
d 



^7 - ^J' 



(10) 



6? = [/-/i + A„ 



which means that the double occupation on the same 
site by the slave particles (holon, spinons and doublons) 
is forbidden. 

In summary, the A/'-component fermion in the reduced 
Hilbert space where the multiple-occupied states are 
truncated is described in slave particle description as fol- 
lows: 



cr-l 



N 



a-1 



cr'=cr+l 

N 



(6) 



Due to the slave-particle constraint (|5|) this representa- 
tion automatically leads to the expected number opera- 
tor: 



cr-l 



N 



^3cr=^]^a^3.cr^Y4.^^'^^.cra'^ Y (7) 



for the holon, spinons and doublons, respectively. \j is 
the Lagrange multiplier for the local constraint. To sim- 
plify the form of i^K, a spinon hopping operator from 
ith to jth site, J-JJ^^, and a creation operator of an an- 
tisymmetric spinon pair between nearest-neighbor sites, 
A^\^'^\ have been defined as 



h- 



(11) 



The spinon pair operator ^^^^^^ is the extension of an 
annihilation operator of a singlet spin configuration for 
the two-component case to a generic A/'-component case. 



III. SINGLE DOUBLON AND HOLON 
PROPAGATOR 

We calculate the single-particle propagator of a holon 
and a doublon based on the Hamiltonian ([8|) and ([9]) in 
the slave-particle representation. 



Atomic limit 



The constraint (|5|) is imposed by a Lagrange multiplier 



Let us start with the atomic limit where J/U = 0. 
Since the atomic Hamiltonian i^at is quadratic in the 
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slave-particle representation, the atomic propagators of 
the slave particles at a jth site are immediately given as 



1 



1 



(12) 
(13) 
(14) 



for a holon, spinons, and doublons, respectively. The 
atomic propagators of spinons and doublons have the 
same form regardless of their species. oOn and denote 
the Matsubara frequency for bosons and for fermions, 
respectively. 



B. Mean-field assumption of spin-incoherence 

We proceed with the finite- but small- J case by making 
a mean-field approximation. In general, due to the effect 
of the hopping Hamiltonian, the system becomes coher- 
ent and exhibits a long-range magnetic order. In order to 
see such phases, the system should reach a temperature 
region lower than the magnetic and charge hopping en- 
ergy scales. However, in the spin- incoherent Mott physics 
case of interest in the present paper, both spin and charge 
coherence are expected to be suppressed due to thermal 
effects. This is a common feature of the atomic limit, 
but to reproduce the finite bandwidth in terms of single 
doublon and holon spectra it is necessary to take into 
account the effect of the kinetic energy Hk- As a simple 
way to describe the spin-incoherent regime, we use the 
assumption, which is valid for J <C ksT <C /7, that 
the spinon propagation is well described by the atomic 
one: 



(15) 



Note that the atomic propagator does not have transla- 
tional symmetry in general. Indeed, Q\^{rj^iuJn) includes 
the local potential coming from the Lagrange multiplier 
Aj, which is potentially site dependent. The mean-field 
treatment of \j is required to recover the translation- 
invariant paramagnetic background in the above frame- 
work. Thus we replace the Lagrange multiplier by the 
homogeneous one: 

A, ^ A. 

Then, the local potentials (p!Q]) also become homogeneous 
by definition: 



e^, e^, and 



e . The mean- 
field A is determined by Eq. (j5j) averaged in the atomic 
limit, 



N{N -I) 



/(e^) + 7V6(6^) = l, (17) 



where f{x) = and b{x) = are, respec- 

tively, the Fermi and Bose distribution functions. Equa- 
tion (pT|) is a saddle-point equation which minimizes the 



atomic-limit free energy. This self-consistent equation 
can be analytically solved for ksT <^ U: 



+ 27V2 - 97V - 6 



-{U-fi)/kBT 



4 



}) 



1/2 



(18) 



For N = 2, it gives back the analytic form given in 
Ref. [41]. As long as J <C ksT <C /7, the atomic limit 
provides the suitable physics. Therefore in the spin- 
incoherent region the above mean-field theory works well 
even if the hopping J is finite. 

Within the mean-field assumption (p!6|) . the atomic 
propagator of the spinons also becomes site independent: 
Gati'^j I'^^n) ^at(^^n)- We cau thus usc the following 
form for the spinon propagator: 



^at(^^n) 



1 



(19) 



At variance with usual mean-field theory, we include here 
the dynamical fiuctuations. The local spin dynamics 
coming from the thermal fiuctuation is thus retained in 
this approximation. 



C. Non-crossing approximation 

Let us consider the full doublon and holon propaga- 
tors, based on the atomic-limit mean-field. To take into 
account i^K at a filling of one or less than one particle 
per site, we use the NCA [s^- This method gives a result 
similar to that from the retraceable path approximation 
by Brinkman and Rice [55|, and is reasonably tractable 
and accurate to describe the physics of single hole mo- 
tion in a MI background. In addition, the NCA allows 
for the control of the chemical potential and tempera- 
ture, which means that one can extend the calculation 
to an inhomogeneous case by the LDA. The NCA dia- 
grams contributing to the self-energy of a doublon and 
holon are shown in Figs. [Hc)-[T]^g). The self-energy for 
(16) doublons is given by two types of diagrams S^ias and 
Hal • Sal comes from the scattering among doublons 
and spinons, and S^lal involves the process in which a 
holon is produced and absorbed. The holon self-energy 
can be constructed in a similar way. The two parts of 
the self-energy diagrams, H^^-*^^ and are illustrated 

in Figs, ma) andlUb). 

In principle, the self-energies are determined by solv- 
ing a set of self-consistent equations for the doublons and 
the holon. The self-energies 

and 

Ecria2 li^k the dou- 
blon and holon propagators as seen in Figs.IIJb) and[Hg). 




FIG. 2. Two examples of diagrams which are not included in 
our NCA: (a) an example of crossing diagrams and (b) one of 
the off-shell diagrams which contains intermediate processes 
at higher energies than the main energy scale, i.e., creation of 
doublons in the holon propagation. The solid, double-solid, 
and dashed lines denote the propagators of the holon, the 
doublon, and the spinon, respectively. Note that the doublon 
and holon propagators shown here mean a bare propagator, 
while all lines denote full propagators in Fig. [H 



FIG. 1. The possible NCA self-energy diagrams for a holon 
[(a) and (b)] and for doublons [(c)-(g)]. The solid, double- 
soUd, and dashed lines denote the full propagators of the dou- 
blon, the holon, and the spinon, respectively. The indices 77 
and T] appearing in the doublon and spinon propagators are 
dummy spin ones, in which the summation over possible spin 
values is taken. The holon self- energies S^(^) and cor- 
respond to (a) and (b). One of the two parts of the doublon 
self-energy Yf^^a-^ includes (c)-(f), and the other Yf^'^^^ is il- 
lustrated by (g). In the simple NCA idea, each propagator 
should be dealt with as the full ones, but in the approximation 
shown in this paper, the spinon propagators are replaced by 
the bare, namely, atomic- limit, ones. In addition, diagrams 
(b) and (g) are off-shell, because they must not be relevant 
in the energy regime ^ U since U is very large. Thus in this 
paper diagrams (a) and (c)-(f) are taken into consideration. 



However, Y^^^^ and T^^}^ can be neglected in the present 
case, as demonstrated below. In the strongly interacting 
case, the two diagrams Fig.[T](b) and (g) are off-shell be- 
cause the intermediate processes creating an additional 
holon and doublon cost an additional energy ~ U . Thus, 
as long as we focus on the physics around the energy scale 
~ [/, the contribution of such diagrams should be negli- 
gible. In particular, this approximation is expected to be 
very good when the system is in a MI state at a filing of 
one particle per site. In Fig. [21 we show examples of the 
diagrams which are neglected in our NCA calculation. 



as 



(20) 



2(7V- 1) 



•77=0-2 + 1 



N 



(72 — 1 



ri=ai + l r)=l 



(21) 



where V is the total number of lattice sites. We have also 
introduced: 



Wb= ^ANzh{€'°)[h{€'°) + l]J, 

= y'8(iV-l)z6(eb)[6(eb) + i]j, 



(22) 



which correspond to the half bandwidths of the lower 
and upper Hubbard bands, respectively, as seen below. 
Because spontaneous breaking of the SU(7V) symmetry 
is unlikely in the spin-incoherent temperature region, the 
single-particle property of the doublons is independent of 
the doublon species: 



Then Eq. (|2T]) is simplified more as 



(23) 



(24) 



Consequently, the self-consistent equations of the self- 
energy of a holon and a doublon are decoupled and given 



which has a form similar to that of Eq. (|2Q|) . The form 
of the self-consistent equations (|2Q|) and (|24|) leads to the 
momentum independence of the self-energies, and thus 
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of the propagators. Using the Dyson equation, the self- 
consistent equation is analyticahy solved, and the result- 
ing propagators are 



Finally, the spectral functions are obtained |56j via ana- 
lytic continuation and are: 



h/d 



4:h 



tkv - eh/d 
H^h/d 



(26) 



The spectral function shows the same semicircular behav- 
ior as for a single hole in a half-filled t-J model discussed 
in Ref. [37] . Note that in that case the slave bosons would 
be condensed as a consequence of the long-range antifer- 
romagnetic order. In addition, this result is similar to 
that of the retraceable path approximation [55| . 

The bandwidth (|22]) as a function of temperature and 
chemical potential is shown in Fig. [3l Wh/Wi depends 
only on N: Wd/Wh = y^2{N - l)/iV, which monotoni- 
cally increases and asymptotically reaches \/2 as N goes 
up. Wd is larger than I^h for N > 2. As in Ref. f41], in 
the SU(2) case, the shape of the two bands is the same. 
Figure [3] shows that the temperature dependence is dif- 
ferent depending on whether or not N = 2. While the 
bandwidth for N > 2 monotonically increases with tem- 
perature, it decreases for = 2 and /Jj/U > 0. 

Let us look at the single-particle properties of the orig- 
inal fermions. Using the slave-particle representation (|6]), 
the Matsubara Green function of the original fermion is 
expressed as 

G.,.'{rj,f,T) = -{Trh\{T)hM){TrhiAr)b],^A^)) 

-X:''^{r.d,,.,(r)4,^,^,(o)) 

77=1 r]' = l 

N N 

77=cr+l r]'—(T'-\-l 

X (T.6],^(r)6,v,,,(0)) 

cr-1 N 

+ E E {T^di.^,{T)d],^^,^{0)) 
X (?V&J,,(r)6,v,^,(0)) 

AT cr'-l 

+ E E<^-^i.''-(^)4.<^v(o)) 

77=cr+l r]' — l 

X (?V&],,(r)6,v,^,(0)), (27) 

where Vj^j' = Vj — Tj/, and denotes the imaginary time 
order. By applying the mean-field assumption (fT9|) and 
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FIG. 3. (Color online) The holon bandwidth Wi,/V^J for 
different chemical potentials as a function of temperature and 
chemical potential: (a) for N — 2^ (b) for A/" = 6, and (c) for 
A = 10, respectively. The lines from upper to lower denote 
IJi/U = 0.5, 0.3, 0.1, 0.0, and -0.1, respectively. 



replacing the atomic propagator in terms of the spinon, 
the propagator of the original fermion, Scr,cr' , t), is 
found to be nonzero only for a = a' and Vj = Vj'. We 
thus take Q^^^f^rjj^^r) = Sjj^dc^^cr'Gai^)- In addition, 
the assumption of SU(A/') symmetry (|23]) simplifies more 
the form of the original fermion propagator. The Fourier 
transform of the propagator is thus written as 



(28) 



and by analytic continuation and Lehmann representa- 
tion, one can obtain the spectral function as 



f{e° -huj) + b{e^) 



Ah{e°/h-uj) 



{N -l)AA{e^/h + Lo). 

(29) 



The spectral functions for A' = 2, 6, and 10 are, respec- 
tively, shown in Figs. HI O and [6l This analytic form 
implies, as mentioned above, that the holon and doublon 
spectra form the lower and upper Hubbard bands, respec- 
tively. Since the centers of the lower and upper bands are 
located, respectively, at = {e^ — e^)/h= {U — /i)/h and 
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FIG. 4. (Color online) The spectral functions of the fermionic atoms for A/" = 2 as a function of the chemical potential for 
different temperatures: (a) ksT < J {ksT /U = 0.025, J/U = 0.05), (b) ksT = J {ksT /U = J/U = 0.05), and (c) ksT > J 
(kBT/U = 0.075, J/U = 0.05). The chemical potential runs from ^/U = to 0.7. 
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FIG. 5. (Color online) The spectral functions of the fermionic atoms for A/" = 6 as a function of the chemical potential for 
different temperatures: (a) ksT < J {kBT/U = 0.025, J/U = 0.05), (b) ksT = J (ksT /U = J/U = 0.05), and (c) ksT > J 
(ksT/U = 0.075, J/U = 0.05). The chemical potential runs from /j^/U = to 0.7. 
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FIG. 6. (Color online) The spectral functions of the fermionic atoms for A/" = 10 as a function of the chemical potential for 
different temperatures: (a) ksT < J {ksT/U = 0.025, J/U = 0.05), (b) ksT = J {ksT/U = J/U = 0.05), and (c) ksT > J 
(ksT/U = 0.075, J/U = 0.05). The chemical potential runs from /j^/U = to 0.7. 



— {e^ — e^)/h = fi/h, the band gap is [/ — (VFh + ^d)- For doublons. 
N > 2, the upper and lower Hubbard bands are always 
asymmetric because of the absence of particle-hole sym- 
metry. In addition, the weight of the doublon band is 
larger than that of the holon, because the possible num- 
ber of doublon states increases with the species of the 
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general it is not easy to directly deal with the effect of 
inhomogeneity, and thus we use the LDA to obtain a 
tractable approximation of (|3T]) . In the LDA framework, 
formula (|3T]) would be identical to the one for the ho- 
mogeneous case (|3Q|) . In what follows, we use Eq. (|3Q|) 
to calculate the DPR spectrum in the same manner as 
discussed in the previous paper [41] in which the inho- 
mogeneity effect of the trap is taken into account by the 
LDA, and the obtained result shows good agreement with 
the experimental data [49|. 



FIG. 7. (Color online) A sketch of optical lattice modulation 
and double occupancy. Due to the modulation perturbation, 
the system is excited, and doubly occupied sites are created. 
In experiments, the number of formed atom pairs is measured 
as a function of the modulation time duration, and the pro- 
duction rate is estimated. 



The DPR is given by the two-particle correlation func- 
tion, which includes vertex corrections. Here we are in 
the strongly interacting regime {J /U <C 1), and we ignore 
the vertex correction as a simple approximation. ^] 



IV. DOUBLON PRODUCTION RATE 

Using the obtained spectral functions ([26|) . we calcu- 
late the DPR spectrum of the optical lattice modulation. 
In the spectroscopy, the amplitude of an optical lattice 
in which the atom cloud is confined is modulated, and 
the created double occupancy is measured as shown in 
Fig. [71 As shown in the Appendix, the DPR per site can 
be obtained from a second-order calculation [53 as 

PD{uj) = -^^uj\mxl{io), (30) 

where XkI*^) is the Fourier transform of the retarded 
correlation function of the kinetic energy Xk(^) ~ 
-ie{t){[TtH^{t),H^{{))]), and 5F is the modulation pa- 
rameter in the lattice model, given by 5F = [dJ/dV — 
dU/dV]SV, where SV is the amplitude of the optical lat- 
tice modulation. 

In order to derive the DPR spectrum formula (|3Q]) . the 
system is assumed to be homogeneous, so the trap is not 
included in the Hamiltonian. It is possible to extend this 
formulation to an inhomogeneous case [47^. Then the 
corresponding response function is replaced by 

(5F)\l{co) ^ -i / dt{[S{t),Sm)- (31) 

The operator S is defined as 5 = {SF)HK — {SU)Hp where 
dU = {dU/dV)SV and is the trap potential term of 
the Hamiltonian. The retarded correlation function is 
computed using the Hamiltonian H + Hp. The above 
formula can be used directly in situations for which a di- 
rect computation of the correlation function in the pres- 
ence of the trap potential can be implemented by use of 
numerics such as Monte Carlo simulations and density- 
matrix renormalization group approaches. However, in 



We now compute the retarded correlation function 
X^(co') for fillings of one or less than one particle per site 
in the spin-incoherent intermediate temperature regime. 
We start with the corresponding imaginary time cor- 
relation function Xk(t) = —{TrHK{r)HK{0)). In the 
same way as in the calculation of the spectral func- 
tion of the original fermions, the analytic continuation of 
the time-ordered correlation function in imaginary time 
leads to the retarded correlation function in real time: 
Xk(^) = XK(^^n cj + iO^), where XK(^^n) is a Fourier 
transform of Xk ('?")• Contrarily to the case of numerical 
evaluations of correlation functions in imaginary time, 
for which there is no straightforward way to perform the 
analytic continuation, here we use our analytic form to 
do so. This is definitely one of the advantages of the 
technique used in the present paper, when computing 
frequency- or time-dependent correlations. 



The result for = 2 is in extremely good agreement 
with the experiment of Ref. [49] , as discussed in Ref . [4l| . 
As we detail below, a similar analytic calculation can be 
also done for the case of A/'-component systems, and this 
allows for a direct comparison to experiments, via the 
LDA for a trapped system. 



The slave-particle representation is useful to clarify the 
physical meaning of the correlation function Xk(t). By 
applying the spin-incoherent assumption ([19]), the corre- 
lation Xk(t) can be written, for fillings of one or less than 
one particle per site, as xkM = XkM + XkM + Xk^M, 
where XkI'^)' XkI'^)' X^i'^) ^^^^ respectively, given 
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(a) 



(b) ^. 






(c) 



(d) 







FIG. 8. Diagrams corresponding to the kinetic-energy cor- 
relation functions: (a) xki^)^ (b) XKi^)^ ^ind (c) Xk^ (''")• The 
sohd, double-sohd, and dashed hnes, respectively, denote the 
propagators of the doublon, the holon, and the spinon. The 
shaded rectangle includes a vertex correction which is not 
considered in this paper. 



as 

X^{t) = NJ^ Yl Qtt{r)gl,{-T){Trh,{T)h]{r)h,{0)^^^^^ 

(32) 

X {TrDl^^^^ {T)Dj,,,a, ir)Dla,,^ (0) A,.4<.3 (0)), 

(33) 



<7l ,(72 ) 

X {Trhl (r)Z)t {T)hj (0) A,,,^, (0)) 



+ (t ^ -t) 



(34) 



They are diagrammatically illustrated in Fig. [8l In or- 
der to simplify the form of the equations the secondary 
doublon operator Dj^^^cr^i which annihilates the doublon 
consisting of a ai- and a <j2-component atom at site j, 
has been introduced: 



(CTI > (72) 
(ai = (72) 
((71 < (72) 



(35) 



The correlation function Xk(t) can be intuitively in- 
terpreted as follows: At initial time, a pair consisting of 
a doublon and a holon is produced by i^K(O), and they 
move in the system. Then the motion of the created 
doublon and holon scrambles up the spin configuration 
of the initial state. For the correlation function to be fi- 
nite in the spin-incoherent case, the spin configuration of 
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FIG. 9. (Color online) A sketch of the contributions to the 
kinetic-energy correlation function. There are four possible 
nearest-neighboring pairs in an equilibrium state: (a) singly 
occupied and empty site, (b) doubly and singly occupied site, 
(c) singly occupied sites and (d) doubly occupied and empty 
site. Based on these configurations of the pair of nearest- 
neighboring sites, we can categorize the three types of contri- 
butions, Xk{^) for ps^'ir {^) ill the left panel, Xk{^) for P^ir 
(b) in the middle panel, and Xk^('^) for pairs (c) and (d) in 
the right panel. The N — 2 case is taken here for the sake of 
simplicity, but it can easily be extended to general N cases. 



the final state must be the same as the initial one. The 
most relevant motion would thus be a retraceable path 
as proposed by Brinkman and Rice [55]. Eventually, the 
doublon and holon go back to the original point of the 
production, and the final state created by acting Hy^{t) 
reproduces the initial state. 

As illustrated in Fig. [9l depending on the spin config- 
uration of the atoms in the initial state, the terms in the 
correlation function, Eqs. ([32|) - ([34|) , contribute as fol- 
lows: XkI'^) nearest-neighboring pairs of singly oc- 
cupied and empty sites, Xy.{^) foi" nearest-neighboring 
pairs of doubly occupied and singly occupied sites, and 
Xk^(t) for nearest-neighboring pairs of singly occupied 
sites and pairs of doubly occupied and empty site. From 
this interpretation, Xk(''') Xk(''') expected to be 
suppressed when the system is at a filling of one particle 
per site. Thus only Xk^('^) leads to important contribu- 
tions to the DPR spectrum. In contrast, in going away 
from the filling of one particle per site, the contributions 
of Xk(^) and Xk(^) appear. 

If one neglects vertex corrections, the two-particle cor- 
relation functions of Eqs. ([32]) - ([34|) are contracted by 
Wick's expansion. As a result, the correlation functions 
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are analytically given as 

xI{t) = ViV^jV(T)C?^-T)C?h(T)gh(-T), (36) 
N{N^ -m + 2,) ^2 



X g''{T)g''{-T)g^{T)g^{-T), 



(37) 



xtir) = -VN{N - l)zJ' 



[g^{-r)]^ g^{T)g\T) 



[g^^rfg'^{-T)g\-T) 



(38) 



By moving from the imaginary-time domain to the real- 
time one by analytic continuation, and by taking the 
imaginary part of the correlation functions, the analytic 
form of the DPR per site at a filling of one or less than 
one particle per site is finally obtained as 



Pd(w) 



{5Ffio [dv 



hu 



27r 



f{hv -hLo)- f{hv) 



^A^{v)A^{v-uj) + 



N{N -\)zP (l + 26(e'')) 
2 

X {f{2e° - hv) - f{2e^ - hv + hio) 



f{hv) + b{2e'') 



X A^{v)A'^{2e°lh-v + oj) 
N{N-l)zJ^ (1 + 26(6*^)) 



X f{2e^ -hv-hu)- f{2e^ - hv) 



f{hv) + b{2e'') 



X A^{v)A'^{2e^/h-v 



(39) 



The DPR spectra for different A^s and chemical poten- 
tials are shown in Figs. [TqHT21 where the small hopping 
is fixed to be J /U = 0.05. To illustrate the temperature 
dependence, the spectra for kBT/U = 0.025 (< J/U)^ 
0.05 (= J/U), and 0.075 (> J/U) are given. As ex- 
pected, the dominant peak is found to appear around 
uj = U/h for every /i//7 and A/", and the peak becomes 
sharper as ja/U gets closer to 1/2, and temperature is 
lowered. Away from filling unity, another small peak in 
the lower frequency regime appears. It occurs because 
Xk becomes relevant due to the hole doping. The spec- 
tral weight of this small peak away from filling unity, 
e.g., at /J./U = 0.1, tends to be suppressed for any N as 
temperature increases. As shown in Fig. [TT] and [121 the 
weight of the peak around uj = U/h for N > 2 increases 
with N. This is due to the enhancement caused by the 
larger spectral weight shown in Figs. |4]-[6]as N increases. 
For such A^s, the spectral weight in the low- frequency 
regime away from /i/U = 0.5, which comes from 



suppressed, in contrast to what happens for the N = 2 
case. However, interestingly, the spectral weight in the 
low- frequency regime for = 6 and 10 also increases 
with temperature, while this strong tendency is not seen 
in the case of A^ = 2. This is due to the finite contri- 
bution of Xk because the doublon band enhanced by the 
larger A^ reaches a; = in such a parameter regime. This 
is shown in Figs. [5] and [6l 

In addition to the plots of Figs. [TQl[T2l we also directly 
fit our results to the experiment [50]. In such experi- 
ments, done with ^^^Yb atoms, the system is expected to 
be dominated by the MI, and thus our calculation scheme 
at a filling of one or less than one particle per site would 
be appHcable, using an LDA calculation to take the trap 
into account. The results, using our theoretical analysis, 
when taking the parameters corresponding to the experi- 
ment are shown in Fig. [131 The experimental parameters, 
namely, the hopping energy J/U, the trap frequency, and 
the modulation amplitude 5F, are taken as follows: (a) 
J/U = 0.0089, 27T X (172, 139,64) [Hz], and SF = 0.085, 

(b) J/U = 0.016, 27T X (175, 141, 69) [Hz], and 5F = 0.09, 

(c) J/U = 0.030, 27r x (181, 146, 78) [Hz] , and 6F = 0.10, 

(d) J/U = 0.061, 27r x (187,151,86) [Hz], and 5F = 
0.115, and (e) J/U = 0.091, 27r x (193, 155, 86) [Hz], and 
6F = 0.125, respectively. The number of atoms in the 
trap is commonly assumed to be 1.87 x 10^. The fol- 
lowing temperatures are determined by the least-square 
fits to the experimental data: (a) ksT/U = 0.0719, (b) 
ksT/U = 0.0577, (c) ksT/U = 0.0909, (d) ksT/U = 
0.0963, and (e) ksT/U = 0.179. Figure [HI shows good 
agreement with the experimental result, which supports 
the validity of our theory. 



V. SUMMARY 

We have computed doublon and holon excitations of 
strongly interacting A^-component fermions in optical lat- 
tices in the spin-incoherent regime. This corresponds to 
a temperature region between the superexchange cou- 
pling and the interaction. As an effective Hamiltonian 
to extract the physics at an energy scale of order ~ /7, 
the symmetric SU(A') Hubbard model has been studied, 
which means that the Hubbard interaction is indepen- 
dent of the internal degree of freedom of the fermions. 
The theory presented in Ref. [4l|pwhich reproduces well 
the experiment with ^^K atoms [49], has been extended 
to an A/'-component fermion case, and the analytic form 
of the single-particle spectral functions for fillings of one 
or less than one particle per site has been obtained. Our 
approach is based on the slave-particle representation 
in which the original fermion operators are represented 
by a fermionic holon, N species of bosonic spinous, and 
N{N — l)/2 species of fermionic doublons. We have em- 
ployed a combination of mean-field theory, a diagram- 
matic approach, and the NCA to take into account the fi- 
nite particle hopping J, and we have captured the physics 
of the hole-doped systems for large interaction J <^U. 
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FIG. 10. (Color online) The DPR spectra per site, PD(uj)/[(6F)^J/h], as a function of modulation frequency for A/" = 2 in a 
cubic lattice {z = 6): (a) ksT/U = 0.025 {ksT < J), (b) ksT/U = 0.05 (/c^T = J), and (c) ksT/U = 0.075 (ksT > J). The 
hopping parameter is taken to be J/U = 0.05. The solid, dashed, and dotted lines denote = 0.5, 0.3, and 0.1, respectively. 
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FIG. 11. (Color online) The DPR spectra per site, Pj:>(uj)/[{SF)^ J/h], as a function of modulation frequency for = 6 in a 
cubic lattice {z = 6): (a) ksT/U = 0.025 {ksT < J), (b) ksT/U = 0.05 (ksT = J), and (c) ksT/U = 0.075 (ksT > J). The 
hopping parameter is taken to be J/U = 0.05. The solid, dashed, and dotted lines denote = 0.5, 0.3, and 0.1, respectively. 
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FIG. 12. (Color online) The DPR spectra per site, Pd{(jo)/[{SF)'^ J/h], as a function of modulation frequency for = 10 in a 
cubic lattice {z = 6): (a) ksT/U = 0.025 (feT < J), (b) ksT/U = 0.05 {ksT = J), and (c) ksT/U = 0.075 {ksT > J). The 
hopping parameter is taken to be J/U = 0.05. The solid, dashed, and dotted lines denote /j/U = 0.5, 0.3, and 0.1, respectively. 



As an application to the calculation of the experimen- 
tal observables, the DPR induced by dynamical periodic 
modulation of optical lattices as a function of modulation 
frequency has been also computed, both for the homoge- 
neous system and for the trapped system, in an LDA. As 
shown in the Appendix, the DPR spectrum as a second- 
order response to the optical lattice modulation is di- 
rectly related to the retarded kinetic-energy correlation 
function. We have discussed the DPR spectrum without 



vertex corrections, and we have presented the analytic 
form constructed by the obtained spectral functions of 
the doublon and the holon. 

From the obtained analytic form in the case of homo- 
geneous systems, we have obtained the DPR spectra as 
a function of temperature, chemical potential, and com- 
ponent number and we have compared the different 
behaviors for the different Ns. While the large peak 
structure around the interaction U exists regardless of 
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the value of some differences have been observed in 
the regime of low modulation frequency. In the compari- 
son, we have focused on two different effects leading to an 
enhancement of the spectral weight in the low-frequency 
regime. The first one is a doping effect: in going away 
from half-filling, the low-frequency spectrum appears as a 
consequence of the system becoming metallic. This effect 
has been found to be suppressed as N increase. The sec- 
ond effect is the temperature: the spectral weight in the 
low- frequency side tends to increase with temperature. 
However, unlike the first effect, we find that the spectral 
weight is enhanced as goes up. Therefore the proper- 
ties of the spectra for different Ns will be most markedly 
different for the low-energy part of the spectrum. 

The theory presented in this paper has several ad- 
vantages: First, the finite-temperature dynamics can be 
dealt with analytically. For such dynamical correlations, 
numerical approaches cannot be straightforwardly ap- 
plied because of the difficulty of numerical implementa- 
tion of the analytic continuation; second, our theoretical 
technique allows for the control of the chemical potential, 
in principle. Note however that our approximations are 
expected to work well at a filling close to the MI state. 
This means that inhomogeneous systems in the presence 
of a trap potential can be also discussed by using an LDA. 
Indeed, in Ref. [41], this approach has been applied to the 
SU(2) symmetric Hubbard model with a harmonic trap 
potential, and quantitatively precise agreement with the 
experiment [49*] has been obtained. 

Using the extension of this approach to trapped sys- 
tems we have compared our results for the DPR spec- 
tra shown in Fig. [131 foi" which the presence of the trap 
potential is taken into account by LDA, with ^^^Yb ex- 
periments. The temperature has been determined by the 
best fit to the experimental data [50], and the obtained 
results for the DPR peak are in good agreement with the 
experiment. 

In recent years, the symmetric SU(A/') systems have 
been being realized in experiments with alkaline-earth- 
metal(-like) ^^Sr @, ^ in addition to ^^^Yb atoms. 
Current fermionic atom systems in such experiments are 
still at high temperature. Therefore our theory is ex- 
pected to work very effectively to compare up-coming 
lattice modulation experiments in such a temperature 
regime. 

Finally, we would like to mention some prospects of 
our study. The first is to apply this technique to the 
calculation of thermodynamic functions such as entropy. 
It is hard to measure temperature directly in experi- 
ments, and the measurement of entropy is used instead. 
Thus by computing the entropy within our theoretical 
framework, we can make a more straightforward com- 
parison with the experiment. The second is to extend 
the theory to general A^-component mixtures away from 
the SU(A^) symmetry limit. Although the SU(A^) sym- 
metry has been assumed throughout this paper, the 
slave-particle representation and the NCA calculation 
would be still applicable away from the SU(7V) symmet- 



ric point. However, the self-consistent equations for the 
self-energies [Eqs. (|2Q|) and (|2T]) ] remain complicated, and 
the issue would be how to solve the self-consistent equa- 
tions. Another prospect is to develop this technique to 
capture the low- temperature physics. The difficulty of 
the application of the present technique to spin-coherent 
systems is that we have here assumed fully incoherent 
spins. Namely, the spinon propagators are replaced by 
the atomic propagators, which means that even nearest- 
neighbor spin correlations are ignored. Thus the key to 
improve the technique for lower temperature would be to 
modify the spin-incoherence assumption (p!5]) . Such an 
improved technique would allow for the crosscheck of the 
theoretical predictions [60|, Ell • 
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Appendix A: Formulation of DPR spectra 

The derivation of the DPR formula is briefly reviewed. 
For simplicity, we only consider the homogeneous case, 
but the more general case including an inhomogeneous 
potential such as a trap can be discussed. Such a gen- 
eral argument can be found in Ref. [47]. We start with 
a generic Hamiltonian of interacting atoms in optical 
lattice potentials defined in D-dimensional continuum 
space, which is written as follows: 

H = Ho^ jdrVo^(r)p{r), (Al) 

where Kp(^) = Yl^=i Vbcos^(/cx^) is the optical lattice 
potential, and Hq is an unperturbed Hamiltonian of in- 
teracting atoms in free space. 

For a deep optical lattice potential {Vq ^ /i), the 
Hamiltonian H is well described by the Hubbard model. 
Then the parameters, the hopping J and on-site interac- 
tion [/, are given as a function of lattice depth Vq. For 
example, if the Wannier function is assumed to be ap- 
proximated as a Gaussian wave function, the hopping J 
and on-site interaction U are estimated [2] as 



J; 



U 



4 ^ (V,V^' 




if) 



3/4 



(A2) 



(A3) 



where and are, respectively, the recoil energy and 
s-wave scattering length of atoms in free space. 



13 



We consider an amplitude modulation perturbation of 
an optical lattice. For deep lattices the modulation ef- 
fect of the lattice potential can be described by replac- 
ing the amplitude of the static lattice potential Vop as 
Vq Vb[l + SV cos{ujt)]. Then the parameters of the 
lattice model also follow the replacement, and the mod- 
ulation parts are derived up to first order in SV: 



U ^U[1^5Ucos{ujt)], 
J ^ J[l + SJcos{ujt)], 



(A4) 
(A5) 



where SU = {d\nU/dVo)Vo and SJ = {d\n J / dVo)Vo . 
In the case of a Gaussian Wannier function ()A2p 
and jM]), 6U ^ 3/4 and 5 J ^ 3/4 - ^/Vo/E^• Thus 
the time-dependent perturbation by lattice modulation 
is written as {dUH\j + S J Hk) cos{ujt) ^ where H\j = 
[/ V • cri>cr2 addltlou, by making use of 

the form of the Hubbard Hamiltonian, the perturbation 
can be rewritten as 5UH cos{ujt) + SFHk cos(co't), where 
6F = 5J — SU. Thus the considered Hamiltonian with 
the lattice modulation is written as 



H{t) = H^ SUH cos{ujt) + SFHk cos{ujt). (A6) 

Extending the doublon number projector for = 2, 
we can define the total number operator of a doublon by 
the Hubbard interaction as 



(A7) 



where we have defined the total doublon number as a 
total sum of all species of doublons. This projection 
operator (|A7p truncates only the empty and singly oc- 
cupied state, and thus it could not count the doublon 
number perfectly for A" > 2 since the projected states in- 
clude multiparticle occupancies of more than three such 
as three- and four- fold occupancy and so on. However, 
because such multiparticle occupations are away from the 
main energy scale in the case considered here , Eq. ()A7p 
would be also identical to the total doublon number in 



the case of multicomponent fermions. The DPR per site 
is defined as a time average of time derivative of A^d over 
a single period of modulation: 



1 



27r/cj 



T+27r/u 



dt 



d (ATp) 
dt V ' 



(A8) 



where (• • •) means the thermodynamic average by the 
Hamiltonian ()A6p . and V is the total number of lattice 
sites. 

We implement second-order perturbation theory in 
term of d{Nj:^)/dt. Then we use the following mathemat- 
ical trick: Using Eq. ()A6p . we rewrite the total doublon 
number operator as 



H{t) -Hk-{H + SFHk) cos(wt) 



(A9) 



From a straightforward calculation up to second order, 
the terms apart from H{t) are found to contribute as os- 
cillatory terms, and they cancel due to the single-period 
time average. Thus Pd('^) can be rewritten as 



27r/w 



L 



(AlO) 



where we have used the identity d{H{t))/dt = {H{t)). 
This is equivalent to the definition of the energy absorp- 
tion rate [46|. This equivalence was numerically estab- 
lished for spin- 1/2 one-dimensional fermions in Ref. (45| . 
The second-order response of the energy absorption rate 
can be calculated by linear response. Therefore one can 
finally obtain the formula as 



^'DH = -^u;ImxgH, 



where Xk('^) ^ 
kinetic-energy-retarded 



2hVU 

Fourier component 



(All) 



of the 

correlation function 
Xlit) = -iO{t){[HK{t),HKm)o where {■ ■ ■)o denotes 
the statistical average by the unperturbed Hamiltonian 
H. 
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FIG. 13. (Color online) The dimensionless DPR spectra 
Pd(^) scaled by J/h as a function of modulation frequency 
huj/U for an = 6 trapped system with different parame- 
ters: (a) J/U = 0.0089, (b) J/U = 0.016, (c) J/U = 0.030, 
(d) J/U = 0.061, and (e) J/U = 0.091. In all cases, 1.87 x 10^ 
trapped atoms and a modulation amplitude SF = 0.08 are 
taken. The solid line and points with an error bar, re- 
spectively, denote the theoretical result and the ^^"^Yb ex- 
periments [50]. The temperatures are determined by the 
least-square fit to the experiment, and the values are (a) 
ksT/U = 0.0719, (b) ksT/U = 0.0577, (c) ksT/U = 0.0909, 
(d) ksT/U = 0.0963, and (e) ksT/U = 0.179. 







